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We present our theoretical and numerical results on thermodynamic properties and the micro- 
scopic mechanism of two successive transitions in vanadium spinel oxides J4V2O4 (A=Zn, Mg, or 
Cd) obtained by Monte Carlo calculations of an effective spin-orbital-lattice model in the strong 
correlation limit. Geometrical frustration in the pyrochlore lattice structure of V cations suppresses 
development of spin and orbital correlations, however, we find that the model exhibits two transi- 
tions at low temperatures. First, a discontinuous transition occurs with an orbital ordering assisted 
by the tetragonal Jahn- Teller distortion. The orbital order reduces the frustration in spin exchange 
interactions, and induces antiferromagnetic correlations in one-dimensional chains lying in the per- 
pendicular planes to the tetragonal distortion. Secondly, at a lower temperature, a three-dimensional 
antiferromagnetic order sets in continuously, which is stabilized by the third-neighbor interaction 
among the one-dimensional antiferromagnetic chains. Thermal fluctuations are crucial to stabilize 
the collinear magnetic state by the order-by-disorder mechanism. The results well reproduce the 
experimental data such as transition temperatures, temperature dependence of the magnetic sus- 
ceptibility, changes of the entropy at the transitions, and the magnetic ordering structure at low 
temperatures. Quantum fluctuation effect is also examined by the linear spin wave theory at zero 
temperature. The staggered moment in the ground state is found to be considerably reduced from 
saturated value, and reasonably agrees with the experimental data. 

PACS numbers: 75.10.Jm, 75.30.Et, 75.50.Ee, 75.30.Ds 



I. INTRODUCTION 

Geometrical frustration in strongly-correlated systems 
is one of the long-standing problems in condensed matter 
physics. Frustration suppresses a formation of a simple- 
minded long-range order and results in nearly-degenerate 
ground-state manifolds of a large number of different 
states. Many well-known examples are found in frus- 
trated antiferromagnetic (AF) spin systems. There, all 
the antiparallel spin conditions between interacting pairs 
cannot be satisfied at the same time because closed loops 
contain an odd number of sites. The degeneracy due to 
the frustration yields nontrivial phenomena such as com- 
plicated ordering structures, spin liquid states, and glassy 
statesi^ Besides the spin degree of freedom, charge or- 
dering phenomena are also much affected by the geo- 
metrical frustration. 3 Quantum and thermal fluctuations 
play important roles in these systems, which are difficult 
to handle in a controllable manner. 

Pyrochlore lattice is a typical example of the 
geometrically-frustrated structures, and it consists of 
a three-dimensional (3D) network of corner-sharing 



FIG. 1: Cubic unit cell of the lattice structure of vanadium 
spinel oxides AV2O4. (a) The pyrochlore lattice of vana- 
dium cations (red balls), (b) The 3D edge-sharing network of 
VO6 octahedra. Oxygen ions are on the corners of octahedra. 
Tetrahedral A sites are omitted. 



tetrahedra as shown in Fig. ^ (a). Spin sys- 

tems on the pyrochlore lattice have been intensively 
studiedi2i4i£i£iL£i2iiii In particular, for quantum S = 1/2 
spin systems with only nearest-neighbor interactions, it is 
predicted that a macroscopic number of singlet states lie 
inside the singlet-triplet gap.- Several types of symmetry 
breakings are predicted within the singlet subspace, e.g., 
dimer/tetramer ordering, but without magnetic long- 
range ordering ^LLS Antiferromagnetic (AF) classical 
spin systems on the pyrochlore lattice are also believed 
to show no long-range ordering at any temperature^^ 
Due to the three dimensionality and the large unit cell 
(16 sites in the cubic unit cell), pyrochlore systems re- 
main a big challenge to theoreticians and are still far 
from comprehensive understanding. 

We can find many pyrochlore systems in real com- 
pounds. Most typically, pyrochlore systems are realized 
in so-called B spinel oxides, where only £?-site cations are 
magnetic in the general chemical formula AB2O4. Fig- 
ure ^(b) shows the B spinel structure, which consists of 
a network of edge-sharing BOq octahedra. B cation is 
shown by ball in center of each octahedron, and octahe- 
dron corners are occupied by oxygen ions, while nonmag- 
netic ^4-site cations are not shown for simplicity. With 
omitting oxygen ions on the corners of the octahedra, one 
obtains the pyrochlore lattice of B cations in Fig. ^ ( a ) • 

In this paper, we will investigate insulating vanadium 
spinel oxides, AV2O4 with divalent A-site cations such 
as Zn, Mg, or Cd. In these compounds, each V 3+ cation 
has two 3d electrons in a high-spin state by the Hund's- 



FIG. 2: Spin ordering structure proposed for ZnV204 on the 
basis of the neutron scattering results. The ordering pattern 
consists of staggered AF chains in the ab plane (red solid lines) 
which stack with a four-times period in the c direction as up- 
up-down-down-- • • (blue dashed lines), (b) Projection of (a) 
from the z direction. Symbols of + and — denote the up and 
down spins, respectively. 



rule coupling, and they are Mott insulators^ Thus, we 
may consider that pyrochlore spin systems with S = 1 
are realized. Curie- Weiss temperature is estimated from 
the magnetic susceptibility as |0cw| ~ 1000K, and any 
long-range ordering does not occur down to significantly 
lower temperatures than |6cw| due to the geometri- 
cal frustration^ For instance, in ZnV204, a structural 
phase transition occurs at T c \ ~ 50K from the high- 
temperature cubic phase to the low-temperature tetrag- 
onal phase with a flattening of VC-6 octahedra in the 
c direction^ Successively, an AF transition occurs at 
T C 2 ~ 40KJ£ Neutron scattering experiments revealed 
that the AF ordering structure below T c2 is a collinear 
one which consists of the staggered AF chains in the ab 
planes stacking in the c direction with a four-times period 
as up- up-down-down-- • ■ as shown in Fig. [2J 13 These 
two successive transitions are commonly seen in com- 
pounds MgV 2 4 and CdV 2 4 ^i& This indicates that 
the degenerate ground-state manifolds in the pyrochlore 
systems are lifted at low temperatures in some manner. 

A few years ago, Yamashita and Ueda proposed a 
scenario to explain the mechanism of the transitions in 
v4V2 04pi£ Their approach is based on a valence-bond- 
solid picture for S — 1 spins and takes account of the 
coupling to Jahn- Teller (JT) lattice distortions. They 
claimed that the first transition at T cl is due to the JT 
effect which lifts the degeneracy of the spin-singlet lo- 
cal ground states at each tetrahedron unit. This sce- 
nario based on the spin-JT coupling is appealing, how- 
ever, some difficulty still remains. The problem is that it 
is difficult to explain the magnetic transition at a lower 
temperature T c2 . In this approach, a finite energy gap 
is assumed between the spin-singlet ground-state sub- 
space and the spin-triplet excitations, and a low-energy 
effective theory is derived to describe a phase transition 
within the spin-singlet subspace. High-energy excitations 
with total spin 5^0 are already traced out at the 
starting point, and therefore their effective model has no 
chance to describe the AF ordering within their theory. 

A similar JT scenario was also examined for classi- 
cal spin systems. 18 In this case, although the problem to 
have AF order does not exist, there is another difficulty 
to explain the following generic difference from chromium 
family of ACr2 04 (A=Zn, Mg, or Cd). These chromium 
oxides are also B spinels and magnetic Cr cations consti- 
tute a pyrochlore lattice. However, in contrast to the two 
transitions in vanadium compounds, the chromium com- 
pounds exhibit only one transition, i.e., the AF order ap- 
pears simultaneously with the structural transition^SS 



This clear difference is generic, being independent of diva- 
lent A cations, and ascribed to the difference of magnetic 
cations V 3+ and Cr 3+ , which cannot be explained by the 
classical spin approach based on the spin-JT effect unless 
there exists essential difference in the model parameters 
between the two families. 

Therefore, these spin-JT type theories appear to be 
insufficient to explain the mechanism of two transi- 
tions in vanadium spinels AV2O4. These insulating 
compounds are undoped states of LiV204 which ex- 
hibits a unique heavy fermion behavior^i^S The ori- 
gin of the mass enhancement is still controversial be- 
tween the scenario based on the Kondo effec1*2 3 *2i2& 
and the scenario of strong correlations with the geo- 
metrical frustrationi2£i2L22i22i31 Since the doping of Li 
shows systematic changes of magnetic^ and transport 
properties^ as well as the phase diagram ; 12 i 32 under- 
standing of undoped materials may give a starting point 
to discuss the doped state in L1V2O4. Therefore, it is 
also highly desired to clarify the mechanism of orderings 
in the undoped compounds AV2O4. 

The generic difference between vanadium and 
chromium spinels mentioned above suggests an impor- 
tance of t 2g orbital degrees of freedom. In the case of 
chromium spinels, each Cr 3+ cation has three electrons in 
threefold t 2g levels and large Hund's-rule coupling leads 
to a high-spin state, and therefore, there is no orbital de- 
gree of freedom. On the contrary, in the case of vanadium 
spinels, since each V 3+ cation has two electrons, the or- 
bital degree of freedom is active. With taking account of 
this t 2g orbital degeneracy, the authors have derived an 
effective spin-orbital-lattice coupled model in the strong 
correlation limit and investigated it by mean-field type 
arguments^ 3 " A reasonable scenario was obtained, but 
discussions were limited to a qualitative level. In order to 
investigate temperature dependences of physical proper- 
ties semiquantitatively accurate enough to be compared 
with experimental data, we need more elaborate analysis. 

In the present study, we will investigate thermo- 
dynamic properties of the effective spin-orbital-lattice 
model derived by the authors in Ref. |33| by extensive 
Monte Carlo (MC) calculations. We will show that this 
model indeed exhibits two successive transitions in a 
reasonable parameter range, and clarify the microscopic 
mechanism of these transitions in detail. First, an orbital 
order appears with the tetragonal JT distortion which 
flattens VC-6 octahedra. This orbital order reduces mag- 
netic frustration partially, and enhances AF spin correla- 
tions in one-dimensional (ID) chains in the ab planes. At 
a lower temperature, the third-neighbor exchange inter- 
action and thermal fluctuations align these ID AF chains 
coherently, and stabilize a 3D collinear AF order. With 
comparing numerical results of physical quantities with 
experimental data, we will show that our theory captures 
essential physics at low temperatures in the vanadium 
spinels AV2O4. 

This paper is organized as follows. In Sec.[Hj we intro- 
duce the effective spin-orbital-lattice coupled model, and 



briefly summarize the mean-field arguments discussed in 
Ref. |3j. Realistic parameter values and MC method are 
also described. In Sec. lIIII we show numerical results in 
comparison with experimental data. We then make sev- 
eral remarks, in particular, on comparisons with experi- 
mental results and other theoretical proposals in Sec. HVI 
Section is devoted to summary. 



FIG. 3: Hopping integrals for the a bonds; t" n and t^ are 
for the nearest-neighbor sites and for the third-neighbor sites, 
respectively. The overlaps between d xy orbitals within the xy 
plane (the solid arrows) and those between d yz orbitals within 
the yz plane (the dashed arrows) are shown. The d zx overlaps 
are similarly taken into account. The gray arrow shows an 
example of second-neighbor pair. 



II. MODEL AND METHOD 

A. Effective spin-orbital-lattice model 

In the present study, we will investigate thermody- 
namic properties of the spin-orbital-lattice coupled model 
which is proposed by the authors in Ref. l83l The Hamil- 
tonian consists of two terms as 



H = H ; 



so 



H 



.IT- 



(1) 



The first term describes exchange interactions in spin and 
orbital degrees of freedom and the second term is for 
orbital-lattice couplings of the Jahn- Teller type. 

The spin-orbital Hamiltonian -ffso is derived from a 
multiorbital Hubbard model with threefold £29 orbital 
degeneracy by the perturbation in the strong correlation 
limits The starting ti g Hubbard model is given in the 
form 

#Hub = J2 E XX^ 7 ** ~ r M a r C jPr + H.C.] 
i,j a.J3 t 

+ 2^ 51 J2 U c*f3, a 'f3>4 aT 4p T ,C % p, T ,C ta , T ,(2) 
i a(3,a.'P' tt' 

where i,j and r, t' are site and spin indices, respectively, 
and a,0 = 1 {d yz ), 2 (d zx ), 3 (d xy ) are orbital indices. 
The first term of 7?Hub is the electron hopping, and the 
second term describes Coulomb interactions, for which 
we use the standard parametrizations^ 

U a p,a'/3' — U 5 aa i5p{l> + Jn(Sap'Spa' + S a pS a 'p'j, (3) 

U = U' + 2J H . (4) 

We do not include here the relativistic spin-orbit coupling 
and the trigonal distortion in the model ©■ Effects of 
these neglected elements will be discussed in Sec. IIVBI 
and E. 

Considering the vanadium spinel oxides are insulators, 
it is reasonable to start from the strong correlation limit 
and treat the hopping term as perturbation. The unper- 
turbed states are atomic eigenstates with two electrons 
on each V cation in a high-spin state. As for the per- 
turbation part, on the basis of the tight- binding fit for 
the band structure^! (see Sec. IIIB|) . we take account of 
hopping integrals of a bonds for nearest-neighbor pairs, 
t™, and for third-neighbor pairs, t^ rd , as shown in Fig. El 
Note that the third-neighbor pair corresponds to the next 
nearest-neighbor pair along each chain. The hopping in- 
tegral between second-neighbor pairs (the gray arrow in 



Fig. |3J) is expected to be small because of the geometry 
of the pyrochlore lattice^ and moreover the exchange 
interaction derived from it is frustrated. 

The second order perturbation in £" n and t^ rd gives the 
Hamiltonian in the form 



H S o = -ffso 
#s n o— 'E 



rr3rd 
W SO > 

dv) 

l o-AF 



dij) 
'o-F 



(i,j) 



rr3rd 
^SO 



dij) 



-J 3 



<(<J» 



'o-AF 



da) 

l o-F 



(5) 
(6) 

(7) 



^^-(A + BSi-Sj) 

X [Tlia(ij)0- - n ja(ij)) + (1 - n ia(ij)) n ja(ij)] , (8) 



"-o-F — ^(1 bi ■ &j)ni a (ij}nj a (ijj, 



(9) 



where Si is the S 



1 spin operator and rij. 



y"L c iaT Ci aT is the density operator for site i and orbital 
a. The summations with (i,j) and ({i,j)) are taken over 
the nearest-neighbor sites and third-neighbor sites, re- 
spectively. Here, a(ij) is the orbital which has a finite 
hopping integral between the sites i and j, for instance, 
a(ij) = 3 (d xy ) for i and j sites in the same xy plane. 
The other parameters in Eqs. ©-© & re determined by 
coupling constants in Eq. J5J) as££ 



J = 


- (tT) 2 /u, 


(10) 


J 3 


= (tl Id ) 2 /u, 


(11) 


A-- 


= (1- V )/(1-3 V ), 


(12) 


B-- 


= 77/(1-377), 


(13) 


C-- 


= (1+^/(1 + 217), 


(14) 


V = 


= Jk/U, 


(15) 



and each site is subject to the local constraint, 
J2 a =i nia ~ 2- Realistic values of these parameters are 
given in Sec. Ill BI 

An important feature of -ffso is the highly anisotropic 
form of the orbital intersite interaction. It is a three-state 
clock type interaction corresponding to three different 
orbital states, in which there is no quantum fluctuation 
since the density operator ni a is a constant of motion. 
This anisotropy comes from the orbital diagonal nature of 
the er-bond hopping integrals which do not mix different 
orbitals. Moreover, the orbital interaction depends on 
both the bond direction and the orbital states in two 
sites. On the other hand, the spin exchange interaction 
is Heisenberg type and isotropic, independent of the bond 
direction. 



FIG. 4: (a) Tetragonal distortion and the level splitting, (b) 
The orbital ordering pattern for the model Q predicted by 
the mean-field argument in Ref.l33l The ferro-type (antiferro- 
type) orbital bonds are shown by the blue solid (red dashed) 
lines. 



The orbital-lattice term Hjt in Eq. QJ reads 



1 ihj) 



(16) 



where 7 is the electron-phonon coupling constant and Qi 
denotes the amplitude of local lattice distortion at site 
i. Here, we take account of only the tetragonal mode in 
the z direction. Note that this simplification breaks the 
cubic symmetry of the system. We choose the sign of Qi 
such that it is positive for a flattening of VC>6 octahedra 
which leads to the level splitting in Fig. 01(a). The second 
term in Eq. (|16[) denotes the local elastic energy of dis- 
tortions. Finally, the third term denotes the interaction 
of JT distortions between nearest-neighbor sites which 
mimics the cooperative aspect of the JT distortion. It 
is reasonable to assume a positive value of A because a 
tetragonal distortion of a VOq octahedron modifies its 
neighboring octahedra in a similar distortion due to the 
edge-sharing 3D network of octahedra in Fig.^(b). For 
simplicity, we here neglect quantum nature of phonons. 
Although JT distortions modify Hso through changes 
of hopping integrals, we neglect these corrections in the 
present study. 

We normalized the variable Qi to absorb the elastic 
constant in the second term in Eq. I|16|l . hence the di- 
mension of Qi is (energy) 1 ' 2 . As a result, the dimensions 
of 7 and A are (energy) 1 / 2 and (energy), respectively. 

For the following discussions, we here briefly summa- 
rize the results of the mean-field type analysis on the 
model Q obtained in Ref. 33. The analysis predicts 
that first, the degeneracy due to the geometrical frustra- 
tion will be partially lifted in the orbital channel. There, 
the anisotropy of the orbital interaction in iJgQ plays a 
crucial role; the interaction is a three-state clock type 
and depends on the orbital states as well as the bond 
direction. The remaining degeneracy of orbital states is 
lifted by the tetragonal JT coupling in Hjt'- A flatten- 
ing of VC>6 octahedra splits the threefold orbital levels 
as shown in Fig. 0] (a), and selects the orbital ordering 
structure as shown in Fig. 0] (b). There, one of the two 
electrons occupies the d xy orbital at every site, and the 
other occupies either d yz or d zx orbital in an alterna- 
tive manner along the z direction. When we consider 
only the nearest-neighbor interactions H^q, this orbital 
occupation induces AF spin interactions on the bonds 
within the ab planes [the blue solid lines in Fig. 0] (b)] 
and ferromagnetic spin interactions on the bonds among 
the ab planes [the red dashed lines in Fig. 0] (b)] . The 



coupling constant for the former AF interaction is JC 
and that for the latter ferromagnetic interaction is —JB. 
Since 77 is a small parameter of the order of 0.1 as es- 
timated in Sec. Ill Bl the former AF interaction is much 
larger than the latter ferromagnetic one. Moreover, the 
ferromagnetic interactions are frustrated for the AF spin 
configuration within the ab planes because of the geom- 
etry of the pyrochlore lattice. Hence, under the orbital 
ordering shown in Fig. |U (b), the AF spin correlations 
develop within the ID chains in the ab planes, and the 
ID AF chains are independent with each other; relative 
angles among the AF moments are not yet determined at 
this stage. The relative angles are partially fixed by in- 
cluding the third- neighbor interactions H^: The third- 
neighbor interactions align the AF moments in two next- 
neighboring ab planes, and lead to a 3D collinear AF or- 
der with the wave vector q = (0, 0, 2n/c). (See Fig. 1101 1 
In the ordered state, however, there are two independent 
AF sublattices; one consists of [110] chains and the other 
consists of [110] chains. The relative angle between the 
AF moments on the two sublattices is still free in this 
mean- field argument. From the spin wave calculation of 
the zero-point energy, we discussed that quantum fluc- 
tuations fix the relative angle and stabilize a collinear 
AF order in Fig. [21 that is consistent with the neutron 
scattering result. 

The mean-field argument gives a reasonable scenario 
for two transitions in AV2O4, but the argument is limited 
to a qualitative level. In order to confirm the scenario and 
understand the experimental results more quantitatively, 
we need more sophisticated analysis, especially for the 
thermodynamic properties of the system. In the present 
study, we will perform the Monte Carlo simulation for 
this purpose. 



B. Parameters 

Here, we estimate realistic values of parameters in the 
model Q which are given by the parameters in the start- 
ing t2 g Hubbard model J2J). As for the hopping parame- 
ters t a p(ri — rj) in i?Hub, a tight-binding fit to the results 
of the first-principle band calculation suggests the dom- 
inant hopping integrals are those of a bonds, and gives 

i™ 0.32eV and i 3rd 0.045eV for nearest-neighbor 

and third- neighbor pair of sites, respectivelyi2L2& For 
Coulomb interactions, there are estimates based on the 
cluster analysis for optical experiments for vanadium per- 
ovskites AVO3, which also have a VOg octahedral unit. 38 
The estimates are U ~ 6eV and Jh ~ 0.68eV, thus, 
i] = J-a/U in Eq. I|15(l is a small parameter of the order of 
0.1. We will set 77 — 0.08 in the following numerical calcu- 
lations. The estimates of t™, i 3rd , and U give J - 200K 
and J 3 ~ 4K, i.e., J 3 /J ~ 0.02. In the Monte Carlo cal- 
culations, we study mainly the case of J3/J = 0.02, but 
we vary the value of J3/J from to 0.05 to examine the 
systematic change by J3. Hereafter, we will set J = 1 as 
an energy unit and the lattice constant of cubic unit cell 



as a length unit (a = b = c = 1), and use the convention 
of the Boltzmann constant k^ = I. 

It is hard to estimate the electron-phonon interaction 
parameters 7 and A, and therefore, we treat them as vari- 
able parameters in the present study. In the present MC 
calculations to confirm the above mean-field scenario, we 
are interested in the parameter region where the orbital 
order in Fig.|3](b) is stabilized by the tetragonal JT dis- 
tortion with a flattening of V0 6 octahedra. The stability 
conditions for this orbital and lattice order will be ob- 
tained in Appendix A by a mean- field type argument. In 
the following MC calculations, we will show MC results 
for typical values of 7 and A which satisfy the conditions 
as 7 2 /J = 0.04 and X/J = 0.15. 



C. Monte Carlo method 

In the present study, we will use MC calculations to in- 
vestigate thermodynamic properties of the effective spin- 
orbital-lattice model QJ. Quantum MC simulations for 
frustrated systems are known to be difficult because of 
the negative sign problem. In the present MC study, we 
neglect quantum fluctuations and approximate the model 
in the classical level. This approximation retains effects 
of thermal fluctuations which may play dominant roles in 
finite-temperature transitions. The quantum nature orig- 
inates only from spin S = 1 operators in the model l|T]l. 
since the orbital interaction is classical and is a diagonal 
one of three-state clock type and since JT distortions are 
also treated as classical variables. Thus, we approximate 
the spin operators by classical vectors with the modulus 
\S\ = 1 [the length of the vector is normalized to give 
the same largest z component S z = 1 (classical part)]. 
Effects of quantum fluctuations will be discussed by us- 
ing the spin wave approximation in Sec. II V CI Thereby, 
the model Q consists of the classical Heisenberg part 
for spins, the three-state clock part for orbitals, and the 
classical phonon part. We use a standard metropolis MC 
algorithm. 

In the actual MC calculations, the MC sampling is per- 
formed to measure spin vectors Si , three-state clock spins 
for the orbital states (defined in Sec. IIII A~2|l . and ampli- 
tudes of the JT distortion Q t at all the lattice sites. We 
typically perform 10 5 MC samplings for measurements 
after 10 5 steps for thermalization. The measurements are 
performed in every ./Vi nt -times MC update, and we typi- 
cally take ATjnt = 2. Results are divided into five bins to 
estimate statistical errors by variance of average values in 
the bins. Here, one MC update consists of several-times 
sweeps (typically twice) for spin directions, orbital states, 
and JT distortions. The one sweep is iV s it e - times trials 
by choosing a site randomly, where N s n e is the number of 
lattice sites. In the sampling on spin directions, we ap- 
ply the so-called pivot rotation when a trial is rejected, 
which is a precession without energy cost: The pivot ro- 
tation of a spin is achieved by a random rotation with 
keeping the relative angle to the mean-field vector de- 



termined by its nearest-neighbor and third-neighbor spin 
and orbital states. This accelerates the MC sampling in 
the configurational space. 

As shown in Sec. Hill the orbital transition accom- 
panied by the JT distortion is first order. To avoid 
a hysteresis and determine the transition temperature 
precisely, we start MC calculations at each tempera- 
ture from a mixed initial condition for orbital and lat- 
tice states in which a half of the system takes a low- 
temperature ordered configuration and the rest takes a 
high-temperature disordered configuration. This tech- 
nique is known to be free from trapping at a metastable 
state for large enough system sizes^S At very low tem- 
peratures, we use an perfectly ordered initial state to ac- 
celerate the convergence. The system sizes in the present 
work are up to L = 12 where L is the linear dimension 
of the system measured in the cubic units, i.e., the total 
number of sites -/V^te is given by L 3 x 16. 



III. RESULTS 

In this section, we present MC results for the model JT} 
in comparison with experimental data. In Sees. IIII At 
C, we show the results for the typical case of J3/J = 
0.02. Systematic changes with J3/J and a generic phase 
diagram will be discussed in Sec. IIII Dl 



A. Two transitions 

In the following, we present MC results for J3/J = 0.02 
to show that the model (JIJ exhibits two phase transi- 
tions with temperature. In Sec. IIII A II we present the 
MC results for the internal energy and the specific heat, 
which show two different anomalies. We also discuss 
the changes of the entropy related to the two transi- 
tions. In Sees. IIII A 21 and 3, we discuss the nature of 
the two transitions by calculating the order parameters. 
In Sec. 1111 A 41 we examine the magnetic ordering struc- 
ture in the low temperature phase, and point out the 
importance of thermal fluctuations. Sec. IIII A~5l contains 
MC results for the uniform magnetic susceptibility for 
comparison with experimental data. 



Internal energy and specific heat 



Figure [S] shows temperature dependences of the inter- 
nal energy and the specific heat per site. The internal 
energy per site is calculated by the thermal average of 
the Hamiltonian (|TJ as 

E=(H)/N sit0 . (17) 

The specific heat is calculated by fluctuations of the in- 
ternal energy as 

(H 2 ) (H)i 



C 



T 2 N si 



(18) 




FIG. 5: (a) The internal energy per site [Eq. 1171 1 and (b) 
the specific heat per site [Eq. 1181 1 at J$/J — 0.02. Error bars 
are smaller than symbol sizes in (a). Typical error bars are 
shown in (b). 



As shown in Fig. 0(a), the internal energy E jumps at 
T ~ 0.19 J. It indicates that a first-order transition oc- 
curs at this temperature. The jump is also found in the 
specific heat at the same temperature in Fig. 03(b). The 
specific heat shows another anomaly at a lower temper- 
ature T ~ 0.115 J. There, we find a systematic enhance- 
ment of the peak as the system size increases, which is a 
sign of second-order phase transition. Thus, MC data in 
Fig. indicate that the system shows two different tran- 
sitions; the first-order transition at To — 0.19J and the 
second-order transition at Tn — 0.115 J. In the following 
sections, the two transitions are to be assigned to the or- 
bital ordering with tetragonal lattice distortion and the 
AF spin ordering, respectively. 

It is noted that the specific heat approaches 3/2 (in 
units of /cb) as T — v as shown in Fig. 0(b). A finite 
value of C at T = is characteristic to classical mod- 
els, and one degree of freedom remaining at the ground 
state contributes 1/2 to C '. Hence, the data in Fig. 
(b) suggests that there remain three degrees of freedom 
at T = in the present model. From the discussions of 
the ordered phase at low temperatures in the following 
sections, they can be ascribed to two transverse modes 
of the AF spin order and one JT mode. 

From the jump of the internal energy in Fig. (a), we 
estimate the entropy jump AS associated with the first- 



order transition from the disordered phase above To to 
the ordered phase below To The two phases have the 
same free energy F = E — TS at To, which gives us the 
entropy difference AS as 



AS = lim 



E{Tp + 5)- E(Tp - S) 
To 



0.4. 



(19) 



This corresponds to ~ 30-40% of In 3 per site. 

The amount of the entropy which is related to the fluc- 
tuation around the second-order transition at Tn is esti- 
mated from the area of the anomalous peak at Tn in the 
plot of C/T as a function of T. Although it is difficult to 
estimate it because of the large system-size dependence 
of C as well as the large error bars in Fig. (51(b), a rough 
estimate is obtained by an interpolation with polynomial 
functions for the normal contribution and a numerical in- 
tegration of the anomalous part. The estimate is roughly 
0.05-0.1, which corresponds to ~ 5-10% of In 3 per site. 

In experiments, the amounts of entropy related to 
the two transitions are also estimated from the specific 
heat. In ZnV2 04, the entropy change in the higher- 
temperature transition at T c \ is ~ 3-4J/mol K, i.e., 
~ 20% of In 3 per V cation4& The entropy related to the 
lower-temperature transition at T c i is small and not es- 
timated quantitatively, but roughly less than U/mol K, 
i.e., ^ 5% of In 3 per V cation4& In MgV204, the former 
is estimated as ~ 3J/mol K, i.e., ~ 16% of In 3, and the 
latter is ~ 0.4J/mol K, i.e., ~ 2% of In 3 per V cation^ 
Our estimates from the MC results in Fig. show semi- 
quantitative agreement with these experimental values, 
and particularly, explain that the entropy related to the 
lower-temperature transition is considerably smaller than 
that for the higher-temperature transition. The small en- 
tropy related to the transition at Tn is likely due to the 
magnetic frustration and a ID AF correlation well devel- 
oped above Tn which will be discussed in Sec. IIIIBI 



2. High-temperature transition: Orbital ordering 

To characterize two transitions in Fig. we calculate 
corresponding order parameters. First, we consider the 
transition at the higher temperature To — 0.19J. Fig- 
ure (a) shows the sublattice orbital moment, which is 
defined in the form 



M 



N si 



E 



iGsublattice 



(20) 



where the summation is taken over the sites within one 
of the four sublattices in Fig. 01(b). Here, Ii is the three- 
state clock vector at the site i which describes three dif- 
ferent orbital states as shown in the inset of Fig. (a) ; 
I L = (1,0) for {xy,yz), I t = (-1/2,^3/2) for (yz,zx), 
and Ii = (—1/2, — V3/2) for (zx,xy) orbital occupations, 
respectively. It is found that the values of Mo for four 
different sublattices have the same value within the error 
bars so that we omit the sublattice index in Eq. I|20|) . As 



shown in Fig.[|)](a), Mo shows a clear jump at the same 
temperature as for the internal energy and the specific 
heat. This suggests that a four-sublattice orbital order- 
ing occurs at To- At low temperatures, Mo approaches 
its maximum value 1, which indicates the four-sublatticc 
orbital order becomes almost perfect there. 

Figures El (b)-(e) show the orbital distribution for four 
sublattices 1-4 shown in Fig. 0] (b) , respectively, which is 
defined as 



iVsit 



E 



(21) 



iGsublatticc 



where a = 1 (d yz ), 2 (d zx ), and 3 (d xy ). The results indi- 
cate that at To the orbital distributions suddenly change 
from equally distributed h a ~ 2/3 in the para phase 
above To to almost polarized n a ~ or 1 for T < To- In 
the orbital ordered phase below To, d yz (a = 1) and d xy 
(a = 3) orbitals are occupied in the sublattices 1 and 4, 
and d zx (a — 2) and d xy (a = 3) orbitals are occupied in 
the sublattices 2 and 3 (Ref. 41). This orbital ordering 
structure is shown in Fig. [7\ This pattern is consistent 
with the mean-field prediction in Fig.|3|(b). 

Accompanying the transition at To, a tetragonal JT 
distortion occurs discontinuously. In Fig. we plot the 
average of the JT distortions which is calculated by 



Q = 52<Qi)/N t 



site- 



(22) 



The positive value of Q below To corresponds to a ferro- 
type tetragonal JT distortion with a flattening of VC>6 
octahedra as mentioned in Sec. Ill Al Therefore, the level 
splitting shown in Fig. 0] (a) is realized. The value of Q 
approaches 2 at low temperatures which is the mean-field 
value obtained in Appendix A. We note that Q is small 
but finite even for T > To- This is because iJjx breaks 
the cubic symmetry of the system as mentioned before. 
Therefore, the discontinuous phase transition at To 
is ascribed to the orbital ordering with the pattern of 
Fig. accompanied by the tetragonal JT distortion with 
the flattening of VC>6 octahedra. From Figs. HH and [S] we 
estimate T Q = (0.19 ± 0.01) J for J 3 /J = 0.02. 
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FIG. 6: (a) The sublattice orbital moment in Eq. 120H at 
J3/J = 0.02. The inset shows the three-state clock vector for 
the orbital state, (b)-(d) Electron density in each orbital for 
four sublattices 1-4 shown in Fig.H(b) for L = 12 [Eq. J2T1 1. 
Error bars are smaller than the symbol sizes. 



FIG. 7: Orbital ordering structure obtained by MC calcula- 
tions. Dark blue and light yellow octahedra in (a) contain V 
cations where (d xy ,d zx ) and (d xy ,d yz ) orbitals are occupied, 
respectively. They stack alternatively in the z direction. Some 
of d zx and d yz orbitals are shown by white and black lobes, 
respectively, (b) is the projection of (a) from the z direction. 
d xy orbitals are singly occupied at all the sites and not shown 



in the figures. 



3. Low-temperature transition: Antiferromagnetic Spin 
Ordering 

Next, we consider the other transition at Tn — 0.1 15 J. 
Figure ED (a) shows the temperature dependence of the 
staggered magnetization defined in the form 

M S = <|/| 2 ) 1/2 , (23) 

where the structure factor is given by 

9 N ^ 

f = J7- E «P(2*«L)E S*(-l) 4y *]. (24) 



This definition looks complicated but nothing but the 
order parameter of the spin ordering pattern shown in 
Fig. [3 Here, the first summation is taken over different 
chains lying in the xy planes (JV c h is the total number of 
the xy chains in the system, i.e., iV c h = 4L 2 where L is 
the linear dimension of the system measured in the cubic 
units), and the second summation ^2 is taken over the 
sites in the i c b-th xy chain, yi is the y coordinates of the 
site i, and If is the z coordinate of the i c h-th xy chain 
measured in the cubic units. We set the normalization in 
Eq. I|24|) such that Mg becomes 1 for the fully saturated 
AF order in Fig. [21 Note that the structure factor / 
cannot be defined only by a real phase factor because of 



Q 




FIG. 8: The average of the JT distortion defined in Eq. J23 
at J3/J = 0.02. Error bars are smaller than the symbol sizes. 



the complicated AF ordering pattern which is expected to 
have the four-period structure in the yz and zx directions 
as shown in Fig-EI See also the discussions in Sec. lIII A4l 
The structure factor is also expressed in a simpler form 
as 



^2gtSi, 



iVsite 

2 

where the form factor gi is given by 

Qi = cos[27r(xi + yi)] + i cos[27r(:Ei - y t )]. 



(25) 



(26) 



Note that gi is specified only by the x and y coordinates of 
the site i because the z coordinate is uniquely determined 
within the cubic unit cell due to the special structure of 
the pyrochlore lattice. [See the projection in Fig. H](b).] 
As shown in Fig. [5] (a), the staggered magnetization 
Ms develops continuously below Tn — 0.1 15 J, and ap- 
proaches the fully saturated value Ms = 1 as T — > 0. 
Figure 0(b) shows the staggered magnetic susceptibility 
obtained by 



K 



Xs 



site 



T 



(d/l 2 )-l(/)l 2 )- 



(27) 



Here, we calculate the susceptibility by using (|/|) in- 
stead of |(/)| in Eq. 1)27(1 for convenience of numerical 
calculations since both quantities agree with each other 
in the thermodynamic limit. The susceptibility xs shows 
a diverging behavior at Tn and the peak value increases 
with the system size. In Fig. 03(c), we show the Binder 
parameter for the staggered magnetization which is de- 
fined as^£ 



,9s 



1 



((1/ 



2\2\ 



3(l/| 5 



(28) 



It is known that the Binder parameter becomes larger 
(smaller) for larger system sizes in the ordered (disor- 
dered) phase, and hence, the crossing point of the Binder 
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FIG. 9: (a) The staggered moment defined in Eq. 12311 . Error 
bars are smaller than the symbol sizes, (b) The staggered 
magnetic susceptibility in Eq. 1)27^ . (c) The Binder parameter 
defined in Eq. 128H . The lines are guides for the eyes. All the 
results are calculated at J3/J = 0.02. 



parameter for different system sizes gives a good estimate 
of the transition temperature. The MC data in Fig.El(c) 
shows the crossing at the same temperature of the di- 
vergence of xs m Fig- ID (b), which indicates the phase 
transition by the order parameter Ms at Tn- 

All these results in Fig. [5] indicate that the phase tran- 
sition at Tn is caused by the continuous growth of the 
staggered magnetic order. From Figs. (b) and (c), we 
estimate T N = (0.115 ± 0.005) J for J 3 /J = 0.02. The 
magnetic ordering pattern will be discussed in the next 
section. 



FIG. 10: Two sublattices; one consists of the [110] chains 
(blue dashed lines) and the other consists of the [110] chains 
(red solid lines). In each sublattice, chains are connected by 
the third-neighbor exchange J3 in the yz and zx directions. 
Black (white) arrows show spins in the sublattice moment Mi 
(M 2 ). 



4- Collinearity of magnetic ordering 

The mean- field argument in Ref. |33J and in Sec. Ill Al 
predicts that the staggered magnetizations develop inde- 
pendently on two sublattices and the relative angle be- 
tween two AF moments is free. The two sublattices are 
shown in Fig.^| one consists of the [110] chains and the 
other consists of the [110] chains. The staggered magne- 
tization Ms defined in Eq. (|23(l is the order parameter 
of the 3D AF spin ordering in Fig. [3 but cannot mea- 
sure the collinearity between two staggered magnetiza- 
tions. For example, Ms takes the value of 1 independent 
of the relative angle between staggered moments in two 
sublattices, Mi and M 2 , shown in Fig. ^3 when both 
sublattice moments saturate. The spin configuration in 
Fig. EH is a noncollinear one, while that in Fig. [2] is a 
collinear one, and both give Ms = 1. 

Here, we measure the collinearity between the stag- 
gered moments in the two sublattices by 



C 12 = ( 



I {M x • M 2 ) 



{M 1 f{M 2 f 



ni) 



(29) 



Here we may consider #12 the angle between M\ and 
Mi- In the present calculations, the two sublattice mo- 
ments are obtained by the real and imaginary parts of 
the structure factor / in Eq. I|24(l as 



Mi = Re(/), M 2 = Im(/), 



(30) 



respectively. If the AF order is collinear, i.e., M\ || M 2 , 
the value of C\ 2 become 1. Figure ^J shows the MC 
results. Below Tn — 0.115J, C\ 2 is rapidly enhanced, 
and becomes larger and approaches 1 as the system 
size increases. This indicates that the magnetic state 
is collinear below T N . Thus, the AF order below T N is 
identified by the collinear AF spin order whose pattern 
is given by Fig. 

The result indicates that thermal fluctuations in the 
classical version of the model Q lift the degeneracy of 
the relative angle between two sublattice magnetizations 
and stabilize the collinear spin structure. This is a kind of 
the so-called order-by-disorder phenomenal In Sec. Ill Al 
and Ref. y3i we discussed that the collinear order may 
appear also by quantum fluctuations. Thus, we can con- 
clude that both thermal and quantum fluctuations favor 
a collinear state of the two sublattice magnetizations. It 
is known that in many frustrated systems, thermal and 
quantum fluctuations favor a collinear state. This is also 
the case for the model Q. 




FIG. 11: The collinearity defined in Eq. |J53J at J 3 /J = 0.02. 
The lines are guides for the eyes. 
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FIG. 12: The uniform magnetic susceptibility in Eq. ()M1|I at 
J3/J = 0.02. Error bars are smaller than the symbol sizes. 



5. Uniform magnetic susceptibility 

We show in Fig. ^| the temperature dependence of the 
uniform magnetic susceptibility calculated by 



iV s , 



X 



T 



-«M t 2 ot ) - <M tot > 2 ), 



(31) 



where M tot is the total magnetic moment per site, M tot = 
I Y2i Si\/N S itc- The result corresponds to the zero- field- 
cool (ZFC) result in experiments rather than the field- 
cool (FC) result since the susceptibility is measured by 
starting the MC simulation from an initial state with 
zero magnetic field at each temperature as described in 
Sec. Ill CI The MC results show a sudden drop at To — 
0.19J and a little continuous change at Tn — 0.115J. 
These features qualitatively agree with the experimental 
results^ 

In experiments, a large difference between ZFC and 
FC results has been found. The difference develops from 
well above T cl , and remains substantial even below T c2 . 
We will comment on this behavior in Sec. IIVDI 
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FIG. 13: The staggered magnetic moment along (a) the xy 
chains [Eq. 1321 1 and (b) the yz chains [Eq. 1331 ] at J3/J = 
0.02. Error bars are smaller than the symbol sizes. 



B. One-dimensional spin correlation in the 
intermediate phase 



As shown in Fig. [§| (a), Ms in finite-size systems are 
enhanced below To even above Tn. In this section, we 
show that there the short-range AF correlation is much 
enhanced along ID chains in the xy planes compared to 
the yz and zx chains. 

The mean-field arguments in Sec. Ill Al and in Ref. [33| 
suggest that the orbital order enhances ID AF corre- 
lations within the xy chains in the intermediate phase 
Tn < T < T> To examine the spatially anisotropic 
spin correlation, we calculate the staggered magnetic mo- 
ments along the chains in different directions. The stag- 
gered moment along the xy chains may be defined by 



M 



(xy) 



1 



E 

Jch=l 



XL 



^'^exp(4 



TTlXi 



1/2 



(32) 

where the summations are taken in the same manner as 
in Eq. 1|24[) . and the phase factor describes the T-i-T-i - ' ' ' 
structure in the xy chains as shown in Fig. [21 Similarly, 
the "staggered" moment along the yz chains is calculated 



M, 



(yz) 



r 2 



AVi, 

E 



-E 

4L^ 



Si exp(27rizi 



1/2 



(33) 

where the first summation is taken over all the yz chains 
in the system and the second summation is taken over 
the sites within the i c h-th yz chain with a phase factor 
of period four in the z direction with considering the "f- 
T-|-|- ■ • ■ structure as shown in Fig. [5] Here, A c h is the 
number of the xy or yz chains, i.e., A c h = AL where 
L is the linear dimension of the system measured in the 
unit cell. The normalization factors in Eqs. 1321) and 1)331) 
are given such that both M g and Mg 
the fully saturated AF state with \Mi\ = 
Eq. EDJ. 

Figure ITSl shows the MC results. The moment in the 
zx chains, Mg , has the same value as Mg Z ' within the 
statistical error bars. Below To — 0.19 J, the moment in 
the xy chains, Mg , is much more enhanced compared 



become 1 in 

|M 2 | = 1 in 



r(v z 



f (zx) 



This 



to that in the yz and zx chains, Mg and M s 
indicates that in the intermediate phase Tn < T < To, 
the AF spin correlations develop mainly within the xy 
chains. 

The system-size dependence of these moments within 
the chains provides further information. Figure IT41 plots 
the system-size extrapolations of the MC data in Fig. 1131 
In the ordered phase below Tn — 0.115J, the MC data 
are extrapolated to finite values, which indicates a long- 
range order. The extrapolated values of Mg and Mg 
to L — ► 00 agree with each other within the error bars as 
expected in this 3D ordered phase. On the other hand, 
in the disordered phase above Tq ~ 0.19J, the MC data 



r( x v) 



r(v z ) 



for both M^ y ' and M§ z ' well scale to \j\fL and are 
extrapolated to zero. (Only the data at T = 0.22J are 
shown in the figures as a typical example but other data 
at T > To show a similar behavior.) This scaling implies 
the exponential decay of the spin-spin correlation along 
the chain as explained below. Compared to the rather 
isotropic behavior for T < Tn and T > To, we find a 



(xy) 



r(vz) 



contrasting behavior between Mg and Mg in the 
intermediate phase Tn < T < Tq; almost all the data of 



(»*) 



M^ yz> well scale to 1/yL even for the smallest system 
size in the present calculations with L = 2, while the 
data of Mg show a clear deviation from the scaling in 
the small-L region. 

In order to analyze this contrasting behavior in the 
intermediate phase, we introduce the correlation length 
£ along the chain, and assume a simple scaling form of the 
spin-spin correlation function; the amplitude of staggered 
spin-spin correlation becomes an almost constant value 
within the length scale £ and decays exponentially for 
further distance than £, that is, 

\{Si'Sj)\ ~ c for r y <£ 

- cexp[-(ry - £)/£] for r l3 > £, (34) 
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FIG. 14: System-size extrapolations of the staggered mag- 
netic moment along (a) the xy chains and (b) the yz chains 
at J3/J = 0.02. Error bars are smaller than the symbol sizes. 
The gray lines are the fits to a linear function of 1/vL. See 
the text for details. 



where c is a constant and r^ = |r^ — Tj\ is the distance 
between sites i and j along the chain measured in the 
cubic units. The staggered moments in Eqs. l|-i2l) and 
(|33[1 can be rewritten by using the spin-spin correlation. 
For instance, Eq. (|32|) is rewritten in the form 



AL 



(xy) 



1 



-1 1/2 



(4L) 
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V2L 






Si ■ Sj cos47r(a;i 

1/2 
(So • S r )\dr 



Vi)) 



(35) 



by introducing the continuum limit to evaluate the sum- 
mation. Note that \J~2L is the length of the chains in the 
L 3 system. Therefore, by assuming the scaling form in 
Eq. QU , we obtain 
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(xy) 
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s/2L 



V2LJ0 

for the case of £ > \/2L, and 
1 r ri 



dr 



1/2 



(36) 
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)} 
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(37) 



for the case of £ < y/2L. The staggered moment in the yz 
chains in Eq. (|33|l is also estimated in a similar manner. 
With considering the above arguments based on the 
scaling form l|34Jl . we can roughly estimate £ from the 
crossover of the scaling from 1/s/Z [Eq. l(3*7jl] to the con- 
stant [Eq. i|36|) ] (+ corrections) in the system-size de- 
pendence in Fig. ^J That is, the linear dimension of 
the system size where the crossover occurs gives a rough 



estimate of £. For instance, 

viates from the scaling of 1/ 

r(y z ) 



at T = 0.13 J, M, 



{xy) 



de- 



L at L £> 8, on the con- 
trary, M^""' obeys the scaling down to L ~ 3. This 
suggests that ^ ~ 6 (~ 16 sites) and ^ ~ 2 (~ 6 
sites). [Note that the longest distance along the chain 
in the periodic L 3 system (iV s it e = 16L 3 ) is L/y/2.] The 
correlation length in the xy direction is about 3 times 
longer than in the yz and zx directions. We note that the 
crossover length systematically shifts to a longer value as 
decreasing temperature, which suggests divergence of £ 
as T — > Tn • [£, is expected to diverge in all the directions 
with the same critical exponent although the divergence 
f £(j/ 2 ) j s no t clear in Fig.^^(b) compared to £( x «).] 

Therefore, we find a ID anisotropy of the spin corre- 
lation: The magnitude of the moment is much enhanced 
and the correlation length is much longer in the xy chains 
than in the yz and zx chains. As predicted in the mean- 
field arguments, the staggered spin correlation develops 
dominantly in the xy chains due to the cooperation of 
the orbital ordering and the geometrical frustration. 

Recently, the neutron scattering experiment has been 
performed in both above and below T c \w^ A clear differ- 
ence of the Q dependence of the inelastic neutron inten- 
sity has been found between the data above and below 
T c i, and ascribed to the ID spin anisotropy due to the or- 
bital ordering. The proposed pattern of the orbital order 
is the same as our result in Sec. IIII A 21 



C. Anisotropy in optical response 

It has been recognized that the interplay between 
orbital and spin degrees of freedom often causes an 
anisotropic electronic state, which is observed in opti- 
cal measurements. For instance, in perovskite vanadium 
oxides, a strong ID nature in the 3D lattice structure has 
been observed in the optical measurement^ and theo- 
retically ascribed to the anisotropic orbital exchange in 
the spin ordered stated Therefore, we expect that the 
ID anisotropy found in the previous section also shows 
up in the optical response. 

We consider the anisotropy of the optical response by 
calculating the spectral weight in the present spinel case. 
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The spectral weight is the total weight of the optical con- 
ductivity defined by 



h = 



2Ti z 



cr MM (w)dw, 



(38) 



where e is the electron charge, fi = h/2it is the Planck's 
constant, ct w is the diagonal clement of the optical con- 
ductivity tensor, and fi = x, y, or z. Here, we consider 
only the optical transfer within the 3d ti g electron bands 
and neglect that to the other bands such as the 3d e g 
bands and the oxygen p bands. In experiments, the inte- 
gral of Eq. Ij38(l is calculated up to an appropriate energy 
cutoff to extract the spectral weight from the relevant 
t2 g levels. The spectral weight in Eq. I|38[) is generally 
given by the kinetic energy in the tig Hubbard model 
Q 47,48 j n ^e strong correlation limit, it is calculated 
by the spin and orbital exchange energy in the effective 
model (0. 46 We show in the following only the results of 
the calculations. Details of the calculations are explained 
in Appendix B. 

In the calculation of the spectral weight for the present 
model, there are three points to be noticed. One is that 
the original Hubbard model (J2J contains not only the 
nearest-neighbor but also the third-neighbor hoppings. 
These two types of hopping contribute differently to the 
spectral weight. The second point is that the directions 
of the electron hoppings are different from the crystal 
axes. Electrons hop along the xy, yz, and zx chains 
which are canted by 45 degree from two crystal axes and 
perpendicular to the rest one. From this, for instance, 
I z is given by the kinetic energy both in the zx and yz 
chains but not in the xy chains. The third point is the 
crystal symmetry. Since the system shows the tetragonal 
lattice distortion below To, there holds a general relation 
in the form 



Iy^h 



(39) 



for T < To . We note that there remains weak anisotropy 
of Eq. 1|39[) even in the high-temperature phase above To 
because of the broken cubic symmetry in Hj^ . With tak- 
ing account of these points, we can derive the following 
expressions 

J * = ^ E [<(^sS)c)+4<(^s 3 o)c)] 5 (40) 



^^zx^xy 



-l 



2K 



site 



E [(( ff sS)c)+4((^Io d )c)] ; (41) 



C=xy,yz 



h = 



-1 



2N f 



silo 



E [<(^S)c)+4<(tf s 3 b d )c)],(42) 



(=yz,zx 



where we set e = h = 1, and (iT?o)c anc ^ (^so')c re P" 
resent the matrix elements of Eqs. © and J7J) in the £ 
chains, respectively (£ = yz, zx, or xy). The derivation 
of Eqs. (gOll-lO will be given in Appendix B. 

In Fig. ED we show the MC results for the spectral 
weight. The spectral weight becomes highly anisotropic 




FIG. 15: The spectral weights for the x, y, and z direction 
calculated by Eqs. |^J-(g5J at J s /J — 0.02. Error bars are 
smaller than the symbol sizes. The solid (dashed) line shows 
T N (To). 



in the tetragonal phase below To with satisfying the re- 
lation (|39|) : The spectral weight in the z direction, I z , is 
suddenly suppressed at To, while those in the xy plane, 
I x and I y , are slightly enhanced there. Moreover, in the 
orbital ordered phase below To, I x and I y increase mono- 
tonically with decreasing temperature, whereas I z shows 
a complicated temperature dependence although the de- 
pendence itself is small; I z slightly decreases above Tn, 
and it turns to increase below Tn- The increase below Tn 
comes from the energy gain in the yz and zx directions 
by the development of the AF spin order. 

The anisotropic electronic state indicated by our MC 
results can be observed in the optical measurements with 
applying the electric field in each direction. It needs a 
clean surface in the single crystal. Unfortunately, it is 
difficult to grow a single crystal large enough thus far, 
and the experimental confirmation of our results remains 
for further study. 



D. Systematic changes for J3/J and phase diagram 

Here, we discuss systematic changes with respect to 
the value of J3/J. Figures ITS1 1171 and 1181 show the spe- 
cific heat in Eq. i|18|) . the uniform magnetic susceptibil- 
ity in Eq. I|31|) . and the orbital and magnetic moments 
in Eqs. (|20|l and (|23|l . respectively, for several values of 
J3/J. The orbital and JT transition temperatures To 
are indicated by the dashed lines in the figures, which are 
monitored by discontinuous changes of C, Xj and Mo as 
well as E and Q (not shown here). The AF transition 
temperatures Tn are indicated by the downward arrows 
in the figures, which are determined by the singular peak 
of C and a continuous development of Ms as well as the 
diverging peak of xs and the crossing point of gs (not 
shown here). 
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FIG. 16: Temperature dependences of the specific heat at 
(a) J 3 /J = 0.0, (b) 0.01, (c) 0.02, (d) 0.03, (e) 0.04, and (f) 
0.05, respectively. The dashed lines (the downward arrows) 
indicate To (7k). Typical error bars are shown. 



FIG. 17: Temperature dependences of the uniform magnetic 
susceptibility at (a) J 3 /J = 0.0, (b) 0.01, (c) 0.02, (d) 0.03, 
(e) 0.04, and (f) 0.05, respectively. The dashed lines (the 
downward arrows) indicate To (Tn). Error bars are smaller 
than the symbol sizes. 



In the case of J3 = 0, the system shows only the or- 
bital and lattice transition, and there is no AF ordering 
down to the lowest temperature. Upon switching on J3, 
the magnetic transition appears at a finite tenrperature 
Tn, below which the 3D AF order exists. As shown in 
Figs. 1161181 Tn increases whereas To slightly decreases 
as J3/J increases. The increase of Tn is easily under- 
stood because the AF spin ordering is stabilized by the 
third- neighbor exchange interaction J3. The reason for 
the decrease of To is not so clear, but it might be due 
to the strong interplay between spin and orbital degrees 
of freedom as well as the frustration between J and J3 
within the same chains. For J3/J <; 0.04, Tn coincides 
with T , and there is only one discontinuous transition 
where both orbital and AF spin moments become finite 
discontinuously. 

Figure 1191 summarizes the phase diagram in the J3-T 
parameter space determined by the analysis of the results 
in Figs. II13II8I as well as other quantities such as E, Q, 
Xs 7 an d 3s ■ The shaded area for Tn < T < To shows the 
orbital ordered phase with the tetragonal JT distortion 
with a flattening of VOg octahedra as shown in Fig.|H(a). 
The ID anisotropy of spin correlation found in Sec. lIIIBl 
appears in this region. The hatched area for T < Tn is 
the AF ordered phase concomitant with the orbital and 
lattice orders. 

From estimates of parameters in Sec. Ill Bl we obtained 
J - 200K and J 3 - 0.02 J - 4K. From Fig. EH these 
estimates give To ~ 40K and T N ~ 20K. The experimen- 
tal values are T cl ~ 50K and T c2 ~ 40K in ZnV 2 4 , 12 
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FIG. 18: Temperature dependences of the orbital sublat- 
tice moment and the staggered magnetic moment at (a) 
J3/J = 0.0, (b) 0.01, (c) 0.02, (d) 0.03, (e) 0.04, and (f) 
0.05, respectively. The dashed lines (the downward arrows) 
indicate To (Tn). Error bars are smaller than the symbol 
sizes. 
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FIG. 19: Phase diagram for the spin-orbital-lattice coupled 
model Q determined by classical Monte Carlo calculations. 
The shaded area Tn < T < To shows the orbital and lattice 
ordered phase. The hatched area below Tn is the AF ordered 
phase concomitant with the orbital and lattice orders. The 
orbital and lattice transition at To (the dashed line) is a first- 
order transition, while the AF transition at Tn (the solid line) 
is a second-order one. The lines are guides for the eyes. 



Thus, two transitions at T c \ and T C 2 in experiments are 
consistently understood by the orbital and lattice order- 
ing transition at To and the AF ordering transition at 
Tn, respectively. The semiquantitative agreement of the 
transition temperatures between our MC results and the 
experimental values is satisfactory with considering the 
assumptions on the derivation of the model Q and on 
parameter estimates in Sec. Ill Bl In particular, we note 
that 7 and A in the JT part in Eq. I|16|) , are parameters in 
our theory for which there has not been any experimental 
estimate to our knowledge as mentioned in Sec. Ill Bl The 
agreement in spite of these assumptions strongly suggests 
that our model (Q captures essential physics in vanadium 
spinel oxides AV2O4. 



IV. DISCUSSIONS 



A. Role of tetragonal JT distortion 

The orbital and lattice orderings at To are considered 
to be caused by the cooperation between the intersite 
orbital interaction in Eq. JSJ and the JT coupling in 
Eq. (|16fl . In this section, we examine the role of the 
JT distortion in more detail. 

We focus on the instability in the orbital sector and 
neglect spins momentarily. That is, we here consider the 
effective orbital model which is derived from the model 
(0 by assuming the spin paramagnetic state in the mean- 
field level as discussed in Ref. 33). We replace Si ■ Sj by 
(Si ■ Sj) = in Eq. JSJ, and obtain the effective orbital 



Ho — Jo 2_^ n ia(i]) n ja(i]), 
(i,3) 



(43) 



where Jo = J (2 A — C). For simplicity, we neglect small 
contributions from the J3 terms in this section. Since 
Jo > 0, the Hamiltonian is an AF three-state clock 
model, which is the same as Eq. (5) in Ref. 33 up to 
irrelevant constants. 

The mean-field argument for the model l|43|) predicts 
that a degeneracy remains partially in the tetrahedron 
uniti^ There are totally 3 4 = 81 different orbital states 
in the four-sites unit, and 30 states among them are in the 
lowest energy state which have four antiferro-type and 
two ferro-type orbital bonds. The 30 degenerate states 
are categorized into two different types shown in Fig. 2 

(a) and (b) in Ref. [33- [The former corresponds to Fig. 0] 

(b) and the latter corresponds to Figs. EH (a) and (b) in 
the present paper.] In the first type, two ferro-type bonds 
do not touch with each other, while they touch at one site 
in the second type. Thus, we expect that in the absence 
of the JT coupling, the orbital degeneracy remains par- 
tially and prevents the emergence of a particular ordered 
state. 

To confirm this prediction, we perform the Monte 
Carlo simulation for the effective orbital model (|43|l . The 
result of the specific heat per site is shown in Fig. 1201 (a) 
(black symbols). The data show a broad peak around 
T ~ 0.4Jo without any singularity or any significant 
system-size dependence. This indicates that the system 
does not show any phase transition, as predicted by the 
above mean- field picture. In Fig. [201(b) (black symbols), 
we plot the temperature dependence of the entropy sum 
defined by 



S(T) 



T C{T) 



T 



dT'. 



(44) 



The entropy sum appears to saturate at ~ 0.5 at high 
temperatures, which is largely suppressed from the ex- 
pected value of In 3 ~ 1.1 for free three-state clock spins. 
This indicates that there remains the entropy at zero 
temperature due to the degeneracy discussed above. 

One way to estimate the zero-temperature entropy So 
is the so-called Pauling's method^ In this method, the 
ground state degeneracy in the whole system is calcu- 
lated simply multiplying the local degeneracy. In the 
present case, the total number of states is 3 Nsit,: , and 
30 of 81 states constitute the degenerate ground-state in 
each tetrahedron as mentioned above. With noting that 
there are N S i tc /2 tetrahedra in the system and assuming 
that the tetrahedra are independent with each other, the 
zero-temperature entropy within the Pauling's approxi- 
mation is obtained as 



So = 



lim 



1 



In 3 
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■In 
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FIG. 20: Temperature dependences of (a) the specific heat 
per site and (b) the entropy sum [Eq. (|44pi ] for the effective 
orbital model Q430 . Black symbols show the results without 
the JT distortion. Gray symbols show the results in the pres- 
ence of the tetragonal JT level splitting with D/ Jo = 0.4 in 
Eq. (gTJ. 



Thus, we obtain the entropy sum as 



1 / 81 
S(T^oc)=ln3-.So = -ln(- 



0.5. 



(46) 



The estimate is close to the saturated value in Fig. |23 
(b). 

The tetragonal JT distortion with a flattening of the 
octahedra as shown in Fig. 0] (a) splits the degenerate 
energy levels, and lowers the energy of one orbital state 
of the first type, i.e., the configuration in Fig.0] (b). In 
this case, all the tetrahedra are in the same orbital states. 
That is, there we expect the orbital ordering and disap- 
pearance of the zero-temperature entropy. The results for 
the case with the level splitting are also shown in Fig. 1201 
(gray symbols). Here, for simplicity, we incorporate the 
level splitting by adding the term 



H 



D 



-E 

3 *-? 



(nil + n l2 - 2n l3 ) 



(47) 



to Eq. (|43[) . which mimics the JT term in Eq. (|16[) . As 
expected, the specific heat in Fig. |20 ( a ) shows a singu- 
larity at T ~ 0.3Jqj which indicates a phase transition. 



Figure l2"01 (b) shows S(T — * oo) ~ In 3, which indicates 
that the phase transition reduces the remaining entropy 
by lifting the degeneracy. 

The results in this section confirm the mean-field pic- 
ture in Ref . '33;, and explicitly show the importance of the 
cooperation between the intersite orbital interaction and 
the JT coupling. 



B. Symmetry of orbital ordered state 

Here, we comment on the spatial symmetry of the or- 
bital ordered state below To- Our result in Fig. [7\ indi- 
cates that the orbital ordering breaks the mirror symme- 
try for the [110] or [110] plane. The symmetry break- 
ing will also appear in the lattice structure through the 
electron-phonon coupling which breaks the d yz and d zx 
symmetry although such coupling is not included in our 
model [the tetragonal JT mode in Eq. (|16fl does not split 
the d yz and d zx levels]. 

Recently, the orbital ordered state with another sym- 
metry has been proposed theoretically. ■ The theory is 
based on the assumption of the dominant role of the rel- 
ativistic spin-orbit coupling. The predicted orbital or- 
der consists of the ferro-type occupation of the d xy and 
{dy Z +idzx) orbitals at every site. This orbital order does 
not break the mirror symmetry. 

The difference of the resultant orbital ordered state 
is important to consider the fundamental physics of the 
present ti g electron system. In ti g electron systems, in 
general, there is keen competition among different con- 
tributions whose energy scales are close to each other; 
the spin and orbital exchange interactions, the JT cou- 
pling, and the relativistic spin-orbit coupling^ In our 
argument, the former two mechanisms play a primary 
role while the last one is not taken into account. On the 
contrary, in Ref. 149 . the relativistic spin-orbit coupling is 
assumed to be dominant. Therefore, the determination 
of the orbital and lattice symmetry is relevant to clarify 
the primary factor in the keen competition. 

Experimental results on the lattice symmetry, however, 
are still controversial. In Refs. 16 and 51, the X-ray scat- 
tering results for polycrystal samples are consistent with 
the symmetry IAi/amd below T c \, which suggests the 
persistence of the mirror symmetry at the low temper- 
ature phase. On the contrary, in recent synchrotron X- 
ray data for a single crystal sample, a new peak is found 
whose intensity is three orders of magnitude weaker than 
a typical main peak4i Unfortunately, the new peak can- 
not conclude whether the mirror symmetry is broken or 
not, but it suggests the different symmetry from the pre- 
vious I4i/amd, i.e., either 74m2 or 14. This controversy, 
probably coming from a small electron-phonon coupling 
in this t<i g system, reveals the necessary of more sophisti- 
cated experiments for larger single crystals. Such experi- 
ments are highly desired to settle the theoretical contro- 
versy. 

Another possible experiment to conclude the symme- 
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try problem is to detect the orbital state directly. This 
may be achieved, for instance, by the resonant X-ray scat- 
tering technique. ■* Besides, another way to distin guis h 
the orbital ordered states in our results and in Ref. |4!j is 
to examine the time reversal symmetry. The (d yz +id zx ) 
orbital order in Ref. |42| breaks the time reversal sym- 
metry even above the magnetic transition temperature. 
A possible experiment is the detection of either circular 
dichroism or birefringence. 



C. Quantum fluctuation effect 

Our classical Monte Carlo calculations have shown that 
with approaching zero temperature, the staggered mo- 
ment increases and saturates at the maximum value, 
Ms — ► 1, when the third-neighbor coupling J3 is finite. In 
physical units, this value should be multiplied by gS/iB 
and corresponds to Ms = 2/ie- Here, /ie is the Bohr 
magneton and the g-factor is assumed to be the standard 
value g = 2. This behavior was plotted in Figs. El (a) and 
1181 (b)-(f). On the other hand, recent experiments show 
Ms ~ 0.6/iB at low temperatures, only less than a half 
of this classical value4i£I It is well known that quantum 
fluctuations due to magnon excitations reduce the am- 
plitude of spontaneous moment in antiferromagnets. We 
expect that this effect is particularly important in frus- 
trated magnets like pyrochlore systems, since frustration 
generally reduces the energy scale of magnon excitations 
and correspondingly enhances quantum fluctuations of 
the staggered moment. 

Here, we examine at T = the effect of quantum fluc- 
tuations of spin degree of freedom on the reduction of 
staggered moment for the model Q, by using the lin- 
ear spin wave theory. For this purpose, we assume the 
perfect orbital order with the ordering pattern in Fig. [7| 
and substitute the density operators in Eq. © by their 
expectation values or 1. That is, we here consider the 
effective quantum spin model on the pyrochlore lattice in 
the form 



Jospin 



= £■* 






7 , ^3^i ' °ji 



(48) 



where the first summation should be taken over nearest- 
neighbor pairs, while the second one be over third- 
neighbor pairs coupled by cr-bond. Under the orbital 
ordering in Fig. [7| exchange coupling constants are fer- 
romagnetic for nearest-neighbor spin pairs on the yz and 
zx chains, J^ = —JB — Jp < 0, and antiferromagnetic 
for those on the xy chains, Jij = JC = Jaf > 0, while 
the third-neighbor couplings are also antiferromagnetic, 
J 3 = J3C. These parameters were defined in Eqs. l(Tn)l- 
(tlTll . It is important that Jaf ^ |</p| 3> J 3 , since fjCl 
and t^ ld <C t 1 ™. Roughly speaking, the small control pa- 
rameter in the spin wave theory is 1/[S x z (number of 
neighbors)]. In the present case, S—l and z=6, leading to 
a control parameter small enough, and therefore we may 
expect quite good estimate from the spin wave theory. 



For the magnetic order determined by the mean field 
arguments in Sec. Ill Al and Ref. 33, we have calculated 
the reduction of staggered moment at T = by using the 
linear spin wave theory. The magnetic structure is shown 
in Fig. [3 and its unit cell contains 8 spins. Details of the 
calculations will be reported elsewhere and here we show 
only the results for the staggered moment. 

When the third-neighbor exchange couplings are zero 
(J3 = 0), magnons have zero energy (zero modes) on 
the planes k x = ±k y in the Brillouin zone. Within the 
linear spin wave theory, magnons are treated as nonin- 
teracting to each other, and these zero modes are excited 
freely without any energy cost. This leads to the loga- 
rithmic divergence of the reduction of moment, AS* — ► 00. 
It is noted that in the pyrochlore spin system in which 
all the nearest-neighbor couplings are antiferromagnetic 
with same amplitude, a half of magnon excitations are 
zero modes throughout the whole Brillouin zone, and this 
results in AS = 00, as a manifestation of strong geomet- 
rical frustration^ 

Once the third-neighbor exchange couplings J3 are 
switched on, zero modes acquire a positive energy and 
the reduction of moment becomes finite. Figurc l2"Tl shows 
the results of the reduction of moment AS as a function 
of J 3 for Jaf = 1 and Jp = —0.113, corresponding to 
i] = 0.08 at which the MC calculations have been per- 
formed. The amplitude of staggered moment is reduced 
to Ms = (S — AS)g^B- In the case of the vanadium 
spinel oxides with S = 1, S - AS > for J 3 ^ 10 -3 Jaf, 
which indicates that the antiferromagnetic order shown 
in Fig. [21 is stable down to very small third-neighbor ex- 
change couplings. For example, at a reasonable value, 
J3 ~ 0.02Jaf, the amplitude of staggered moment is 
Ms ~ (1 — Q.5)g/j,B ~ 0.5g/x B . With using g = 2, this 
corresponds to Ms ~ l^B- Hence, the reduction is sig- 
nificantly large in the present spin-orbital-lattice coupled 
system, and particularly the staggered moment rapidly 
decreases as J3 in the realistic small- J3 region. We con- 
sider that the estimate of Ms is satisfactory and can 
explain the experimental value by carefully tuning the 
model parameters. 

We also note that there may be a small contribution 
from quantum fluctuations in the orbital degree of free- 
dom. In our spin-orbital- lattice coupled model p|l. this 
contribution is not taken into account since we consider 
only a bond for hopping integrals. Other hopping inte- 
grals lead to small orbital exchange interactions, causing 
small quantum fluctuations of orbitals. This will slightly 
reduce the orbital sublattice moment and may lead to a 
further reduction of the staggered spin moment. 



D. Difference of ZFC and FC susceptibility 

We have shown that the MC results for the uniform 
magnetic susceptibility in Sec. IIII A 51 qualitatively ex- 
plain the experimental data in the zero-field-cool (ZFC) 
measurement. In experiments, there has been reported 
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FIG. 21: Reduction of the staggered moment due to quantum 
fluctuations (open symbols) estimated by the linear spin wave 
calculation at Jf/Jaf = —0.113. The staggered moment is 
also plotted (filled symbols), which is given by Ms = (S — 
AS)g/j,B and by assuming g = 2. The lines are guides for the 
eyes. 



that the ZFC and field-cool (FC) data show a large differ- 
ence at low temperatures^ The difference between the 
ZFC and FC data emerges at a higher temperature than 
the structural transition temperature T c \, and develops 
on cooling. One of the mysteries is that the difference re- 
mains finite and large even in the AF spin ordered phase 
below T C 2'. The FC value becomes nearly twice of the 
ZFC one at the lowest temperature. This strongly sug- 
gests that the difference cannot be explained as a simple 
spin-glass phenomenon. 

A key observation is that the starting temperature of 
the difference between the ZFC and FC data as well as 
the magnitude of the difference appears to depend on 
samples. In Ref. Ha. the difference emerges at T ~ 100K, 
and becomes ~ 0.8 x 10 3 emu/V-mol at T cl . On the other 
hand, in Ref. |5l], the difference starts at T ~ 70K, and 
becomes ~ 0.3 x 10 3 emu/V-mol at T cl . We also note 
that the temperature dependences of the ZFC data show 
a sample dependence, especially near two transition tem- 
peratures. These sample dependences imply the impor- 
tance of quenched disorder. Effects of disorder in the 
present spin-orbital-lattice coupled system are interest- 
ing and open problems. 

In our MC calculations, it is difficult to obtain the FC 
results because the system shows a discontinuous tran- 
sition at To: There appears a huge hysteresis when we 
change the temperature successively by using the MC 
sample in the previous run as the initial state in the 
next run. On cooling, the system remains in the high- 
temperature para state well below the true transition 
temperature To &s a metastable state because of an ex- 
ponentially large energy barrier in the first-order tran- 
sition. One way to avoid this hysteresis is to apply 
more sophisticated MC method such as the multicanon- 
ical technique^ This interesting problem is left for fur- 
ther study. 



E. A-site substitution effect 

We have shown that our spin-orbital-lattice coupled 
model Q exhibits two transitions which well agree with 
experimental results in AY2O4. Here, we discuss a dif- 
ference among the compounds with different A cations 
observed in experiments. 

Among j4V 2 4 with A=Zn, Mg, or Cd, CdV 2 4 shows 
a different behavior near T c \ from others^ As temper- 
ature decreases, the magnetic susceptibility for the com- 
pounds with A=Zn and Mg shows a sudden drop at T c \ 
as in our numerical result in Fig. 1121 whereas that for 
the Cd compound shows a sharp increase. Moreover, the 
transition temperature T c \ is substantially higher in the 
Cd case than in the Zn and Mg cases; T c \ = 97K for Cd 
while T c i = 52K and 64. 5K for Zn and Mg, respectively. 
On the other hand, the AF magnetic transition temper- 
ature T C 2 is almost the same among three compounds; 
T c2 = 44, 45, and 35K for Zn, Mg, and Cd, respectively^ 

Our results show that the transition temperature To, 
which corresponds to T c i, is determined by the en- 
ergy balance between the intersite orbital interaction in 
Eq. (JjjJ and the elastic energy gain in Eq. (|16|) . The A- 
site cations are nonmagnetic and locate at the tetrahedral 
sites which are separated from the pyrochlore network of 
V cations, and the change of the A cation may affect the 
lattice structure since AO4 tetrahedra share oxygen ions 
with VOe octahedra. Thereby, we may consider that the 
Cd cation, which has relatively large ionic radius, may 
modify the lattice structure and have an influence on the 
JT energy gain as well as the form and magnitude of the 
orbital interaction. 

Actually, although the symmetry is the same for the 
three compounds, the so-called u parameter of the spinel 
structure is considerably larger for CdV204 than for Zn 
and Mg compounds. The u parameter represents the 
magnitude of the trigonal distortion. The absence of the 
trigonal distortion gives u — 0.375 by definition, and a 
larger u denotes a larger trigonal distortion. For the Zn 
and Mg cases, u = 0.385 and 0.386, respectively, whereas 
u becomes 0.394 for the Cd compound: The trigonal dis- 
tortion is substantially larger in the Cd case than in the 
Zn and Mg cases. This suggests that the peculiar behav- 
ior in CdV204 may be ascribed to effects of the trigonal 
distortion. 

The agreement between the experimental data in Zn 
and Mg compounds and our numerical results, which do 
not include effects of the trigonal distortion, indicates 
that the small trigonal distortion does not play an im- 
portant role in these two compounds. To understand the 
behavior of the Cd compound, i.e., the sharp increase in 
the magnetic susceptibility as well as the higher T c i, it is 
interesting to examine the effect of the trigonal distortion 
in our theoretical framework. In addition to the contri- 
bution from the electron-phonon coupling to the trigo- 
nal distortion, we may have to include more complicated 
hopping integrals not only the cr-bond type but also, for 
instance, 7r-bond type or the second-neighbor hoppings, 
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since the trigonal distortion affects the network of VC>6 
octahedra: These additional hoppings modify the spin 
and orbital intersite interactions in the derived effective 
model in Eq. (JSJ . The complicated effects of this trigonal 
distortion are left for further study. 



V. SUMMARY 

In the present work, we have investigated the micro- 
scopic mechanism of two transitions in vanadium spinel 
oxides AV2O4 with nonmagnetic divalent A cations such 
as Zn, Mg, and Cd. We have focused on the role of 
the tig orbital degree of freedom as well as spin in these 
strongly correlated electron systems on the geometrically 
frustrated lattice structure, i.e., the pyrochlore lattice 
consists of the magnetic V cations. We have derived the 
effective spin-orbital-lattice coupled model in the strong 
correlation limit from the multiorbital Hubbard model 
with explicitly taking account of the ti g orbital degen- 
eracy. The effective model describes the interplay be- 
tween orbital and spin, and reveals the contrasting form 
of the intersite interactions in two degrees of freedom; 
the Heisenberg type for the spin part and the three-state 
clock type for the orbital part. The anisotropy in the or- 
bital interaction originates from the dominant role of the 
cr-bond hopping integrals in the edge-sharing configura- 
tion of VC>6 octahedra. The Jahn- Teller coupling with 
the tetragonal lattice distortion is also included. Ther- 
modynamic properties of the effective model have been 
investigated by the Monte Carlo simulation, which is a 
classical one to avoid the negative sign problem due to 
the geometrical frustration. Quantum corrections are ex- 
amined by the spin wave approximation. Main results are 
summarized in the following. 

Our effective spin-orbital-lattice coupled model ex- 
hibits two transitions at low temperatures. As tempera- 
ture decreases, first, the orbital ordering transition occurs 
assisted by the tetragonal Jahn- Teller distortion. This 
is a first-order transition. Successively at a lower tem- 
perature, the antiferromagnetic spin order sets in. This 
transition is second order. 

For realistic parameter values, our numerical results 
agree with experimental data semiquantitatively. The es- 
timates of two transition temperatures are To ~ 40K and 
Tn ~ 20K, while the experimental values are T c \ ~ 50K 
and T C 2 ~ 40K. The changes of the entropy at the tran- 
sitions are comparable to the experimental values: In 
particular, the small change at Tn compared to the con- 
siderable change at To is well reproduced. The mag- 
netic ordering structure in the low temperature phase 
is completely consistent with the neutron scattering re- 
sults. The magnitude of the staggered moment at T = 
is largely reduced by quantum fluctuations, and the es- 
timate by the spin wave theory reasonably agrees with 
the experimental data. The temperature dependence of 
the uniform magnetic susceptibility is similar to the ex- 
perimental data. The agreement strongly indicates that 



our effective model captures the essential physics of the 
vanadium spinel compounds AV2O4. 

Our results give an understanding of the mechanism 
of the two transitions in the vanadium spinel oxides. 
The present numerical study has confirmed the mean- 
field scenario in our previous publication^ and more- 
over, given more detailed and quantitative information. 
The first transition with orbital and lattice orderings is 
induced by the intersite orbital interaction which is three- 
state clock type and spatially anisotropic depending on 
both the bond direction and the orbital states in two 
sites. The anisotropy lifts the degeneracy due to the ge- 
ometrical frustration inherent to the pyrochlore lattice. 
The tetragonal Jahn- Teller coupling also plays an im- 
portant role to stabilize the particular orbital ordering 
pattern. The obtained orbital ordering structure is the 
antiferro type which consists of the alternative stacking of 
two ferro-type ab planes; (d xy , d yz ) orbitals are occupied 
in one plane, and (d xy ,d zx ) orbitals are occupied in the 
other. Once the orbital ordering takes place, the orbital 
state affects the spin exchange interactions through the 
spin-orbital interplay and reduces the magnetic frustra- 
tion partially. As a consequence, the antiferromagnetic 
spin correlation develops mainly in the one-dimensional 
chains in the ab planes. We found typically about three- 
times longer correlation length in the xy direction than in 
the yz and zx directions. At the second transition, these 
one-dimensional chains are ordered by the third-neighbor 
exchange interaction to form the three-dimensional an- 
tiferromagnetic spin order. It was found that thermal 
fluctuations stabilize the collinear state by the order-by- 
disorder type mechanism. 

There still remain several open problems on this topic 
as discussed in Sec. II" VI the controversy on the symme- 
try of the orbital ordered state, the large difference of the 
zero-field-cool and field-cool data of the magnetic suscep- 
tibility, and the A-site dependence, which is probably re- 
lated to the trigonal distortion. They need further inves- 
tigations from both experimental and theoretical view- 
points. Besides them, another important problem is the 
carrier doping effect on the Mott insulating materials, 
such as the Li doping in Zni_ a; Li a ;V204. In experiments, 
it is known that the doping rapidly destroys the ordered 
states at x = 0, and replaces them by a glassy stated 
Finally, at x = 1 , the system becomes metallic and shows 
a heavy-fermion like behaviori2i*22i We believe that our 
present results give a good starting point to study the 
interesting doping effects. 
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Appendix A 



FIG. 22: Orbital ordering patterns (a), (b), (c), and (d) for 
the categories 1, 2, 3, and 4 specified by Eqs. 15511 . 1561 . 1571 . 
and 15811 . respectively. The ferro-type (antiferro-type) orbital 
bonds are shown by the blue solid (red dashed) lines. 



In this Appendix, we examine the conditions for 7 and 
A in the Hamiltonian (|16fl to stabilize the orbital order in 
Fig. 21 (b) accompanied by the tetragonal JT distortion 
with a flattening of VOg octahedra as shown in Fig. 0] 
(a) at low temperatures. We employ the mean-field type 
argument to discuss the instability at high temperatures 
as in Ref. 33. Here, we consider only the nearest-neighbor 
interactions and neglect small contributions from -fffo ■ 
With considering a spin disordered state and replacing 
Si ■ Sj by (Si ■ Sj) = in Eq. QJ, we obtain the effective 
Hamiltonian for the orbital and JT parts as 



#0-JT — Jo2_^ nia(i]) n 3a{tj) 



H. 



JT j 



(49) 



where Jo = J (2 A — C). The first term is equivalent 
to Eq. IJ4"3|) and to Eq. (5) in Ref. |3J up to irrelevant 
constants. 

We consider the orbital state and the JT distortion in 
the tetrahedron unit shown in Fig. |H(b). For an orbital 
configuration at the four sites, the expectation value of 
the JT part -ffjT for one tetrahedron is written in the 
following quadratic form 



ejT(Q) = Q ■ MQ + Q-A + A-Q, 



(50) 



where Q is a four-dimensional vector which denotes the 
amplitude of the JT distortion in each site as Q =* 
{Qi, Q2, Q3, Q4), and M is the 4x4 matrix whose matrix 
elements are given by Mij = 1/2 for i = j and Mij = — A 
for i ^ j. The vector A describes the electron-phonon 
coupling part and depends on the orbital state; for in- 
stance, A =* (-7/2, -7/2, -7/2, -7/2) for the orbital 
occupation in Fig. |U(b). Here, we consider the physical 
situation in which the system does not show any sponta- 
neous lattice distortion without the coupling to electrons 
of the second and third terms in Eq. J5UJ). This is sat- 
isfied when the matrix M is positive definite. Three of 
the eigenvalues of At are (1 + 2A)/2 and the last one is 
(1 — 6A)/2, and therefore, by noting that we consider only 
positive A as mentioned in Sec. Ill Al the condition is 



0<A<i 
6 



(51) 



The quadratic form of ejx in Eq. H5U|I can be trans- 
formed to the following form 



e JT (P) = (P + B) ■ (P + B) - B ■ B, 



(52) 



where P = LQ, B = L~ X A, and L = M 1 / 2 is well- 
defined since M is a positive-definite matrix under the 



condition of Eq. (|51|l . Therefore, we find that the JT 
energy takes its minimum value 



mm 
e JT 



B ■ B = A ■ M~ l A, (53) 

for the amplitude of the distortions 

gmin = _M~l A _ (54) 



For instance, for the configuration in Fig.0|(b), we obtain 
e™ n - -2 7 7(1 - 6A) for Q min =* (Q,Q,Q,Q) with 
Q = 7 /(1-6A). 

All the different orbital configurations with (n, Q ) = 
or 1 are classified by the energy value ejx into five groups, 
which are described by the vector A as 

Ai = '(-7/2,-7/2,-7/2,7), (55) 

A 2 = '(-7/2,-7/2,7,7), (56) 

A 3 = '(-7/2,7,7,7), (57) 

A 4 = '(7,7,7,7), (58) 

^5 = '(-7/2,-7/2,-7/2,-7/2), (59) 

where arbitrary permutations of components in each vec- 
tor give the same energy. Typical orbital configurations 
are shown in Fig. [221 (The configuration for A5 is shown 
in Fig. Kb).) 

On the basis of the above consideration on the JT en- 
ergy, we calculate the expectation value of the Hamil- 
tonian (|49|) . for the five types of orbital configurations. 
With considering the orbital interaction energy from the 
first term in Eq. I|49l) , the minimized values of the energy 
per tetrahedron are obtained as 



pi 



P-2 



C3 



P-i 



P.5 



4Jo- 
4Jo- 
6J0- 
8J0- 
4Jo- 



7 2 (7-40A) 
2(1 + 2A)(1-6A)' 

7 2 (5 - 26 A) 
(1 + 2A)(1-6A)' 

7 2 (13-28A) 
2(1 + 2A)(1-6A)' 

8 7 2 
1-6A' 

27 2 
1-6A' 



(60) 
(61) 
(62) 
(63) 
(64) 



Note that there are twelve bonds per one tetrahedron. 

By comparing e 5 with other four values e, (i = 1 — 4), 
we obtain the stability conditions for the orbital config- 
uration of Fig. 01(b). First, the conditions e$ < e\ and 
es < 62 are satisfied when 



1 . 1 

To <A< 6- 



(65) 
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FIG. 23: The parameter region of Eqs. I|65|l - H67|| is shown 
by the shaded area. The dashed and solid lines denote the 
conditions of Eqs. 1661 and 1671 , respectively. 



In this range, the condition es < e 3 leads to 
4J (1 + 2A)(1-6A) 



7 2 < 



(66) 



9(1 -4A) 

and the condition e§ < e^ leads to 

7 2 <2J (l-6A)/3. (67) 

The shaded area in Fig. [23 shows the parameter re- 
gion where all the conditions of Eqs. (|65(l - (|67|l are sat- 
isfied. We performed several MC runs to confirm the 
present mean-field type analysis. In the MC calculations 
in Sec. |ffll we take j 2 /J = 0.04 and X/J = 0.15 as a 
typical set in this region. 



Since the pyrochlore lattice consists of three different 
chains in the xy, yz, and zx directions, and we here con- 
sider only the cr-bond hopping integrals which are finite 
only along the chains, K(k) is given by the summation 
of the contributions from each chain as 

K(fc)=^[ Mc ^"(fc)+« c ^ d (fc)], (70) 

C 

where £ = 1 (yz), 2 (zx), 3 (xy) denote the three dif- 
ferent types of the chains, and u^ , v^ are the nearest- 
neighbor and the third-neighbor hopping integrals along 
the C chain, respectively. In the Hamiltonian l|2|l. 1/1 = 
u 2 = u 3 — t™ and v\ = v 2 = v 3 = £^ rd , but they are de- 
liberately denoted by different parameters for later use. 
Each term in Eq. I|7()|l is defined by 

n™(k) = 2cos(fc • <5oi)[4(fc)c u (fc) + c t 11 (fc)c i(fc)] 
+2cos(fc^ 23 )[c t 21 (fc)c 31 (fc)-|-c^(fc)c 21 (fe)], (71) 

K f d (k) = 2cos(2fc • «5oi)[4(fc)c i(fc) + c t 11 (fe)c 11 (fc)] 
+2cos(2fe-<J 23 )[c 21 (fe)c 2 i(fe) + 4 1 (fe)c3i(fe)] I (72) 

and so on, where d a b = da — fib- 
Let us now evaluate the derivatives of K(k) in Eq. JSSJ. 
By using Eq. (|70|l . it is straightforward to obtain, for 
instance, the derivative in terms of k x in the form 



d 2 K(k) 
dkl 



+ 



1 
16 

1 



u 2 K™(k)+u 3 n™(k)] 
[v 2 4* d (k)+v 3K f d (k)]. 



(73) 



Appendix B 

In this appendix, we derive the expressions for the 
spectral weights in Eqs. 140|l - (|42|l for our effective model 

CD- 

The spectral weight, which is the total weight of the 
optical conductivity, is given by the kinetic energy of the 
system (/-sum rule) 42*^ We may safely extend the for- 
mula for the single-band Hubbard model; 



/„ 



2 

7T./n 



i(u)du> 



N nx 



E 
k 



d 2 K(k)\ 



dkl 



), (68) 



to our multiorbital case. (We set e = % = 1, and iV u . c . 
is the number of unit cells.) Here, K(k) describes the 
kinetic energy term of the multiorbital Hubbard Hamil- 
tonian as 

J2 K(k) - E E E t aa Mk)cL(k)c bf 3(k), (69) 



k 



k a > b a >0 



which is the Fourier transform of the first term of Eq. J2J>. 
Here, a, b = 0, 1,2,3 are the atom indices in the primi- 
tive unit-cell containing four V atoms, whose relative po- 
sitions are labelled by 6 = (0,0,0), d 1 = (0,1/4,1/4), 
S 2 = (1/4, 0, 1/4), S 3 = (1/4, 1/4, 0), respectively. 



The derivatives in terms of k y and k z are obtained in a 
similar manner. Thus, with noting N uc — N s - lte /4:, the 
spectral weights for the multiorbital Hubbard model @ 
are obtained in the forms 



4 = — 



l 



4K 



site 



{((HT) ZX ) + ((H™) x% 



+4{{(H$? d ) zx ) + ((#£<%) 



Iu = — 



4K 



site 



{<(#K n W + <(#K n ) 



y~ 



4N S : 



4{<(^ru + ((^rw 



{((#£" W + ((#£")** 



+4{(( H k)v*) + (( H k)**) 



(74) 



(75) 



(76) 



where (H™)^ and (J?k )c are the nearest-neighbor and 
the third-neighbor matrix elements of the kinetic term in 
Eq. J2J along the £ chain, respectively. 

In the strong correlation limit, the kinetic energy cor- 
responds to the spin and orbital exchange energy as dis- 
cussed below. The expectation values of the kinetic terms 
in Eqs. (|74|I - H76|I can be represented as coupling-constant 
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derivatives; 



«*™ = («C^Hu b ) = ^i<ffHu b >,(77) 



(m 



3rd\ 



i) 



d\a.\v^\ 



(H 



Hub/- 



(78) 



Here, we have used the Hellman-Feynman theorem to ex- 
change the order of derivative and expectation value. In 
the strong correlation limit, we can replace the deriva- 
tives in terms of uq and v^ of (-ffHub) by the derivatives 
in terms of the exchange J of the expectation values of 
the effective spin-orbital Hamiltonian (J5J as 







d\n\uc\ 



(Hnub) — 2 







51n|(J) ( 



(TJso) = 2<(i/ s n S) c )(79) 



d 

d In \v^\ 



(Hn u b) 



d 

! ain|(J 3 ) C | 



(H, 



so; 



2<(#f r o d )$0) 



where (J)q = uf/U and {Jz)c, = vf/U are the nearest- 
neighbor and the third-neighbor exchange coupling con- 
stants along the (, direction, respectively; and (-ffgo)c an d 
(^so d )c are the nearest- neighbor and the third-neighbor 
matrix elements of Eqs. (JSJ and J7J) along the £ chain, 
respectively. In the model (J5J, (J)^ = J = (t° n )/U and 
(^s)c = ^3 = (^o- rd )/^7 ancl they are independent of the 
direction £. By substituting all the expectation values 
in Eqs. lfft |l -l|75 |l by Eqs. (JTSJ), (f5Uj>. and the similar ex- 
pressions, we obtain the formulas of the spectral weights 
for the effective spin-orbital coupled model in the strong 
correlation limit as given in Eqs. 1401) - (|42|) . 
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